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This paper deals with maximum entropy completion of partially specified block-circulant ma- 
trices. Since positive definite symmetric circulants happen to be covariance matrices of stationary 
periodic processes, in particular of stationary reciprocal processes, this problem has applications in 
signal processing, in particular to image modeling. Maximum entropy completion is strictly related to 
maximum likelihood estimation subject to certain conditional independence constraints. The maximum 
entropy completion problem for block-circulant matrices is a nonlinear problem which has recently been 
^ solved by the authors, although leaving open the problem of an efficient computation of the solution. 

in 

\q The main contribution of this paper is to provide an efficient algorithm for computing the solution. 



(N 



Simulation shows that our iterative scheme outperforms various existing approaches, especially for 



[- — large dimensional problems. A necessary and sufficient condition for the existence of a positive definite 

^_l circulant completion for unitary bandwidth and block-size is also provided. 



k> I. Introduction 

H 

We consider the problem of completing a partially specified block-circulant matrix under 

the constraint that the completed matrix should be positive definite and block-circulant with an 
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inverse of banded structure. As shown in [4 J a block-circulant completion problem of this kind 
is a crucial tool for the identification of a class of reciprocal processes. These processes, ll22l . 
E51l . Il27ll . are a generalization of Markov processes which are particularly useful for modeling 
random signals which live in a finite region of time or of the space line (think for example 
of an image). In this paper we consider stationary reciprocal processes for which we refer 
the reader to [|24|. |[T3l and references therein. In particular, stationary reciprocal processes of 
the autoregressive type can be described by linear models involving a banded block-circulant 
concentration matrixjj whose blocks are the (matrix-valued) parameter of the model. 

This problem fits in the general framework of covariance extension problems introduced by A. 
P. Dempster [11 J and studied by many authors (see |f!2|, 11201, [[Toll, Il28l IfTSl, lfH, |[T7I, Ifl4l and 
references therein). A key discovery by Dempster is that the inverse of the maximum entropy 
completion has zeros in the positions corresponding to the unspecified entries in the original 
partially given covariance matrix, a property which, from now on, will be referred to as the 
Dempster property. For the special case of block-Toeplitz matrices the inverse of the maximum 
entropy extension turns out to have a banded structure. See |[T2| and ifTTll for a general formulation 
of matrix extension problems in terms of banded algebras and for a thorough discussion of the so- 
called band-extension problem for block-Toeplitz matrices. It turns out that the block-Toeplitz 
band extension problem can be solved by factorization techniques and is essentially a linear 
problem. This is unfortunately no longer true when a block- circulant structure is imposed (6) 
to the extension. 

A relevant fact is that, even when the constraint of a circulant structure is imposed, the inverse 
of the maximum entropy completion maintains the Dempster property. In fact, it has been first 
noticed in Hi for a particular data structure and then proved in complete generality in 0, that if 
the data are compatible with a block-circulant structure then the maximum entropy completion 
turns out to be automatically block-circulant. Otherwise stated, the solution of the Maximum 
Entropy block-Circulant Extension Problem (CMEP) and of the Dempster Maximum Entropy 
Extension Problem (DMEP) with data consistent with a block-circulant structure, coincide. This 
fact permits to simplify the formulation of the entropy optimization problem considerably. Note 
that this property does not hold, for example, for arbitrary missing elements in a Toeplitz structure 

'i.e. the inverse covariance matrix, also known as the precision matrix. 
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since in such a case the completion is not automatically Toeplitz unless the given data lie on 
consecutive bands centered along the main diagonal lfT2ll . 

Since the solutions of the CMEP and of the DMEP with circulant-compatible data coincide, all 
methods available in the literature for the DMEP can, in principle, be employed to compute the 
solution of CMEP In particular, it has been shown that if the graph associated with the specified 
entries is chordal ( lfP9l ). the solution of the DMEP can be expressed in closed form in terms of 
the principal minors of the sample covariance matrix; see 0], [fT5l , 11261 . The sparsity pattern 
associated with the given entries in our problem is not chordal. For non-chordal graphs the 
maximum entropy completion has to be computed iteratively. A straightforward application of 
standard optimization algorithms is too expensive for large size problems, and several specialized 
algorithms have been proposed in the literature; see lUTI . 11281 . fl29), 11231 . These algorithms, 
however, deal with the general unstructured setting of Dempster. In the present work, we propose 
a modified matricial gradient descent algorithm for the solution of the CMEP which naturally 
follows from the variational analysis in 01 and exploits in an essential way the circulant structure 
of our problem. This algorithm compares very favourably with the algorithms proposed in the 
literature for the solution of the DMEP. 

A further issue addressed in the present paper is the feasibility of the CMEP. In O, a 
description of the set of all positive definite completions of a partially specified block-circulant 
matrix has been provided, while in [4] a sufficient condition on the data for a positive definite 
block-circulant completion to exist has been derived. Here we provide a necessary and sufficient 
condition for the feasibility of the CMEP for unitary bandwidth and block-size. 

Necessary and sufficient conditions for the existence of a positive definite completion have 
also been discussed by Grone, Johnson, Sa and Wolkowicz [20] . Although our results may seem 



different they are in fact not in contrast with those in [20J. See Remark 6.2 in Section VI 



The outline of the paper is as follows. In Section |TT] we introduce some notation and state the 



entropy maximization problem. Section III contains a brief review of some of the most popular 



methods for the solution of the DMEP. In Section IV the matricial gradient descent algorithm 
is introduced. An experimental comparison between the proposed algorithm and the methods 
available in the literature is presented in Section |VJ Finally, in Section VI a necessary and 



sufficient condition for the feasibility of the CMEP is discussed. A numerical example which 



sheds light on the feasibility issue is presented in Section VII 
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II. Notation and preliminaries 

All random variables in this paper, denoted by boldface characters, have zero mean and finite 
second order moments. It is shown in 01 that a wide-sense stationary M m -valued process y is 
stationary reciprocal on [0, N] if and only if its covariance matrix, say S^v, has a block-circulant 
symmetric structure, i.e. S^ is of the form 



J iV 



S S x 

Si s 



E T 



s7 

Si 



Si 

s r 



s T 



s7 

So 



s7 



•■■ s7 

Si Sq 



where S fc = Ey(t + A;)y(t) T . We refer the reader to [|9]| for an introduction to circulants; an 
extension of some relevant results for the block-case can be found, for example, in Q. Here 
we just recall that the class of circulants is closed under sum, product, inverse and transpose. 
Moreover, all circulants commute and are simultaneously diagonalized by the Fourier matrix, 



see (21 ) below, of suitable size. 



Let (5 at denote the vector space of real symmetric matrices with N x N square blocks of 
dimension m x m. Moreover, let X b be the set of pairs of indices consistent with a banded- 
symmetric block-circulant structure of bandwidth n, i.e. the set of the (i, j)'s such that 

if | i — j? | < ran =3- (i,j) G X h 

if (i,j) el b ^ (j,i) el b 

if (i,j) Gli^ ((j + m) mod mN , (i + m) mod mN J G I b . 



(an example of this structure is shown in Figure la). The differential entropy H(p) of a probability 
density function p on MJ 1 is defined by 



H{p) = — \og(p(x))p(x)dx. 

JM." 



(2) 
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In case of a zero-mean Gaussian distribution p with covariance matrix H N , it results 

Hip) = \ log(det Sjv) + ^n (1 + log(2vr)) . (3) 

The Maximum Entropy covariance extension problem (CMEP) for block-circulant matrices can 
be stated as follows. 

max {log det Sat | Sjv G & N , *£ N > 0} (4a) 

subject to : 

eiT, N ej = nj, for (i,j) el b and r„ G 7£ 6 (4b) 

Sjy is block-circulant (4c) 



where e^ is the row-vector e^ 



and TZb denotes the set of the given 



... 1 ... 
data consistent with a banded- symmetric block-circulant structure, i.e. 

'R-b '■= {rij G K | (i,j) G 2&, ry = r^ry = T(j+m) modmN (i+ m ) raodmN } • 

We are interested in the case where the ry's represent the entries of a partially specified block- 
circulant covariance matrix, say Rjv- If we remove constraint ( |4"c] ), we get the maximum entropy 
covariance selection problem studied by A. P. Dempster (DMEP). 

Notice that, although in Problem [4] we are maximizing the entropy functional over zero-mean 
Gaussian densities, we are not actually restricting ourselves to the case of Gaussian distributions. 
Indeed, the Gaussian distribution with (zero mean and) covariance matrix solving ([4]) maximizes 
the entropy functional ([2]) over the larger family of (zero mean) probability densities whose 



covariance matrix satisfies the boundary conditions (j4bj), ( |4c| ), see [4, Theorem 7.2]. 

The first question we need to address is feasibility of the CMEP. In [4], a sufficient condition 
for the existence of a positive definite block-circulant completion has been provided, while a 
characterization of the set of all positive definite completions of a partially specified block- 
circulant matrix has been derived in Q. In the present paper, we provide a necessary and 
sufficient condition for the feasibility of the CMEP for unitary bandwidth and block-size. For 



the sake of continuity in exposition, we postpone this issue until Section VI 



III. Algorithms for the covariance selection problem 

Before discussing some of the algorithms proposed in the literature for the covariance selection 
problem, a brief digression into some basics of graph theory is needed. First of all, it is natural to 
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describe the pattern of the specified entries of an Nm x Nm partial symmetric matrix M = (m^) 
by an undirected graph of Nm vertices which has an edge joining vertex i and vertex j if and 
only if the entry m^ is specified. Since the diagonal entries are all assumed to be specified, we 
ignore loops at the vertices. The undirected graph will be denoted by Q = (V,E), where V is 
the vertex set and E the edge set which consists of unordered pairs of distinct vertices. In any 
undirected graph we say that 2 vertices u, v E V are adjacent if (w, v) E E. For any vertex set 
S C V, consider the edge set E(S) C E given by 

E(S) :={(u,v) EE\u,vES} 

The graph G(S) = (S,E(S)) is called subgraph of Q induced by S. An induced subgraph G(S) 
is complete if the vertices in S are pairwise adjacent in Q. In this case, we say that S is complete 
in Q. 

Definition 3.1: A clique is a complete subgraph that is not contained within another complete 
subgraph. 

The sequence vo,Vi, . . . ,Vk with (fj,fj+i) E E, Vi ^ Vj for i ^ j is said a simple path (of 
length k) in Q. Similarly, the sequence Vo, v±, . . . , Vk, Vq is called a simple cycle in Q. A chord 
of a path (cycle) is any edge joining two nonconsecutive vertices of the path (cycle). 

Definition 3.2: An undirected graph is chordal if every cycle of length greater than three has 
a chord. 

Finally, we define the complementary graph of Q = (V, E) as the graph Q with vertex set 
V and edge set E with the property that (u, v) E E if and only if u ^ v and (w, v) <£ E. As 
anticipated in the Introduction, if the graph of the specified entries is chordal ([19]), the maximum 
determinant matrix completion problem admits a closed form solution in terms of the principal 
minors of the sample covariance matrix (see [[3, [fT5l , Il26ll ). However, the graph associated 
with a banded circulant sparsity patterns is not chordal, as it is apparent from the example of 



Figure lb Therefore we have to resort to iterative algorithms. For the applications we have in 
mind, we are dealing with vector-valued processes possibly defined on a quite large interval. A 
straightforward application of standard optimization algorithms is too expensive for problems of 
such a size, and several specialized algorithms have been proposed in the literature ([[TT|. [|28l . 
||29l , [|23l ) which deal with the general, unstructured setting of the DMEP. In his early work 
( ifTTlO . Dempster proposed two iterative algorithms which however are very demanding from a 
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Fig. 1 : Banded pattern representing the given entries for the CMEP in a block-circulant structure 
with N = 12, n = 3, m = 1 (on the left) and associated graph (on the right). The graph is not 
chordal since, for example, the cycle {1,4,7, 10} does not have a chord. 

computational point of view. Two popular methods are those proposed by T P. Speed and H. T 
Kiiveri in Il28ll . that we now briefly discuss. 

a) First algorithm: As mentioned in the Introduction, for the class of problems studied 
by Dempster, the inverse of the unique completion which maximizes the entropy functional has 
the property to be zero in the complementary positions of those fixed in Y, N . Thus, a rather 
natural procedure to compute the solution of the covariance selection problem seems to be the 
following: iterate maintaing the elements of Sjy in Z& at the desired value (i.e. equal to the 
corresponding elements in the sample covariance matrix) while forcing the elements of E^ 1 in 
XI to zero. This procedure can be formally described by Algorithm [T] 
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Algorithm 1 First algorithm (Speed and Kiiveri [28 1) 



Compute all the cliques c t in the complementary graph Q, say {c t , t = 1, . . . , n^}; 

Initialize H N = Hn', 



while some stopping criterion is satisfyied do 
for all the cliques c t in Q do 



end for 
end while 












where 



J Af 



is the mN x miV zero matrix which equals 

-i 



diag 



i(*-l)N-l 

'TV J 



i(*-l)N-l 

J N ) 



ct 



in the positions corresponding to the current clique c t (given a mN x mN matrix M and a set 
a C [1, . . . , Nm), M a denotes the submatrix with entries {rriij : i, j G a}). Every cycle consists 
of as many steps as the cliques in the complementary graph Q (the graph associated to the 
elements in T£). At each step, only the elements in Sjy corresponding to the current clique c t 
(i.e. only a subset of the entries in Z£) are modified in such a way to set the elements of S^ 1 
in the corresponding positions to the desired zero-value. Throughout the iterations, the elements 



,(*) 



in Sy are fixed over Zj,, while the elements of I T, N 



(t) 



vary over Z£. The crucial step in the 
algorithm involves going from S^ - to Y, N and relies on the following Lemma, whose proof 
can be found in ([28, Lemma 2, (i)]). 

Lemma 3.1 (Speed and Kiiveri): Let Q, R and B be positive definite matrices. Then for 

a C V the matrix 



Q- 1 = R- 1 H- 

is positive definite and satisfies 

Q{a, 0) = B{a, 0) 
Q- 1 (a,P) = BT x {a,P) 



,-i 



(B a y l - (R a 



-1 











if a E a and /3 E a 
if a ^ a and (3 ^ a 



(5) 

(6a) 
(6b) 
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b) Second algorithm: The role of T, N and S^ 1 can also be swapped, yielding an alternative 
procedure, which is the analog of iterative proportional scaling (IPS) for contingency tables 11211 . 
Let <f ( H N J be the mN x mN zero matrix which equals 



((^^-((sroj 



in the positions corresponding to the current clique c t in Q (the graph associated with the given 
entries). The second algorithm reads as follows. 

Algorithm 2 Second algorithm (Speed and Kiiveri Il28l0 
Compute all the cliques c t in Q, say {q, t = 1, . . . , n Ct }; 

Initialize T, N = In™', 

while some stopping criterion is satisfyied do 
for all the cliques c t in Q do 

(e«)- 1 = (e<5- 1, )" 1 + p (eS- ,> 

end for 
end while 

Every cycle consists of as many steps as the cliques in the graph of the specified entries Q. 
At each step, only the elements in E^ 1 corresponding to the current clique c t (i.e. only a subset 
of the entries in Z&) are modified in such a way to set the elements of Ejv in the corresponding 



positions to the desired value, namely equal to the sample covariance R^ (see (|6a|)). Through 
the iterations the elements in ( S^ j are fixed over Z£ while the elements of £)/ vary over 

Algorithms [T] and [2] can be seen as special cases of a general procedure whose convergence is 
proved in ll28l . An intuitive, geometric interpretation is also possible. To see this, let us denote 
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with C (C) the class of the cliques of Q (Q) and let us define the "subspaces" 
V L;C = {P > | P c = L c } 

Qm,c = {Q > | (Q _1 )c agrees with M 5 except on the diagonal} 
V L , c = f]V LtC 

c&C 



q m ,c = n q ^ 



c<=C 

where, for every cC [1, . . . , mN], P c is the submatrix P c := {P(i,j) \ i E c, j E c}. In the limit, 
S^ G Vr n ,c H Qj q. The first algorithm begins with Y, N ' G V-r, n ,c and cycles throught c EC 



I-projecting the current estimate of H N onto V-R N ,c^Qi mN ,c while the second begins with £ 



(0) 

N 



Qj £ and cycles throught c EC X-projecting ([7]) the current estimate onto Vr NiC fl Q 1 ^. 
(The fact that the operation performed at each step is an X-projection follows from [28, Lemma 
1 part (ii)]). It follows that, in a suitable sense, both algorithms are an application of the von 
Neumann's alternating projections theorem for computing the projection onto the intersection of 
nonorthogonal subspaces. 

c) Comparison of the two algorithms: The choice of which algorithm is to be preferred 
depends on the application and is very much dependent on the number and size of the cliques in 
Q and Q. In our setting, the complexity of the graph associated with the given entries depends on 
the bandwidth n. In particular, for bandwidth not too large with respect to the completion size 
(which is the case we are interested in) the complexity of the graph associated with the given 
data Q is far lower than the complexity of its complementary (which, for small n, is almost 
complete), see Figure [3] The execution time of the two algorithms has been compared for a 
completion size N = 30 and a bandwidth n varying between 2 and 8. The results are shown in 
Figure [2] and Table |I} It turns out that for n small the second algorithm (which, from now on, 
will be referred to as IPS) runs faster than the first, and thus has to be preferred. 

d) Covariance selection via chorda! embedding: In [8], Dahl, Vanderberghe and Roy- 
chowdhury propose a new technique to improve the efficiency of the Newton's method for the 
covariance selection problem based on chordal embedding: the given sparsity pattern is embedded 
in a chordal one for which they provide efficient techniques for computing the gradient and the 
Hessian. The complexity of the method is dominated by the cost of forming and solving a system 
of linear equations in which the number of unknowns depends on the number of nonzero entries 
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added in the chordal embedding. For a circulant sparsity pattern, it is easy to check that the 
number of nonzero elements added in the chordal embedding is quite large. Hence, their method 
does not seem to be effective for our problem. A preconditioned conjugate gradient method has 
also been proposed in (flU), but the comparison is not carried out in the present paper. 



SK1 vs SK2 - Execution time 




N (completion size) 



Fig. 2: Comparison between the execution time of the first and second algorithm for iV = 30, 
m — 1, n — {1, . . . , 8}. 



IV. A Matricial Gradient Descent Algorithm 
Let XJn denote the block-circulant shift matrix with N x N blocks, 

I m ... 



u 



jV 



L 





i m o 





-'m 
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(a) Q for n = 2 



(b) S for n = 2 





(c) 5 for n = 5 



(d) 5 for n = 5 





(e) Q for n = , 



(f) 5 for n = 



Fig. 3: Graph Q associated with the given data (on the right) and its complementary Q (on the 
left) for iV = 20 and bandwidth n = 2, 5, 8. 
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First algorithm 


Second algorithm 


n 


# of cl. (max. cl. size) 


CPU time [s] 


# of cl. (max. cl. size) 


CPU time [s] 


2 


4608(10) 


9.7877 


30(3) 


0.4109 


3 


2406(7) 


4.1515 


30(4) 


0.1783 


4 


1241(6) 


1.9419 


30(5) 


0.3153 


5 


706(5) 


1.0525 


30(6) 


0.5535 


6 


445(4) 


0.6258 


30(7) 


0.9854 


7 


295(3) 


0.4145 


30(8) 


1.7477 


8 


175(3) 


0.2480 


30(9) 


3.0665 



TABLE I: Execution time of the first and second algorithm for N — 30, m — 1, bandwidth 
n = {2,...,8}. 



let T n G & n +i denote the Toeplitz matrix of boundary data 



Si S S x 



and let E n be the iV x (n + 1) block matrix 

~I m 

U 1 777 



E n 











... l n 



(7) 
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so that, specializing the formulation in Q, the maximum entropy band extension problem for 
block-circulant matrices can be stated as 

max {log det Sat | Y, n G & n , T, n > 0} (8a) 

subject to : 

E n Sat -En = T n , (8b) 



\J N 'E N \J N — H N . (8c) 



where we have exploited the characterization of block-circulant matrices as invariants under the 
similarity, namely C N = U^CatUat. 

Problem ([8]) is a convex optimization problem since we are minimizing a strictly convex function 
on the intersection of a convex cone (minus the zero matrix) with a linear manifold. We shall 
briefly recall the approach to this problem by duality theory following |@]|. To this aim, consider 

the linear map 

A : G n+1 xG N -»■ G N 



(9) 

(A, e) m. E n AEZ + v N eul-e 



and define the set 



£+ :={(A,6) G (6 n+1 x G N ) | (A, 6) G (ker(^)) ± , 

(E n kEl + V N Q\J T N - 0) > 0}. (10) 

C + is an open, convex subset of (ker(A)) ± . Letting (A, B) := trAB T , the Lagrangian function 
of this problem can be written 

L(V N , A, 9) := - tr log Sjv + < A, (EjB N E n - T n ) > , + <0, (U^U^ - £*) > 

= -trlogS w + tr(£ n A£ n T S ]V ) - tr (AT n ) + tr (U^U^) -tr(eEjv) 

and its first variation (at 11 N in direction <5£at G &n) is 

«JL(Ejv, A, 9; SH N ) = - tr (E^Eat) + tr (£ n A£^ Sjv ) 

+ tr((u A reu A ;-e)5s J v). 

Thus SL(E N , A, 9; <$£jv) = 0, V^ G 6jv if and only if 

S^ 1 = E n KEl + UA.9UT - 9. 

July 14, 2011 DRAFT 



DRAFT 



15 



It follows that, for each fixed pair (A, 0) 6 £+, the unique *E° N minimizing the Lagrangian over 

G N ,+ := {Sat 6 &n, £jv > 0} is 

i:° N = {e u aeJ + v N ev T N -ey 1 . 

Moreover, computing the Lagrangian at S^r = £yv results in 

L(-E° N , A, 6) = -trlog ([E n XE T n + U^GU^ - 0) 



(ID 



-l 



tr 



{E n KEl + UnQVI - 0) 



-i 



tr(AT n ) 



(E n AEj + u N eu T N -ey 

= trlog (E n KEl + U^BU^ - 0) 
+ \xI mN - tr (AT n ) . 

This is a strictly concave function on C + whose maximization is the dual problem of (CMEP). 
We can equivalently consider the convex problem 

min{J(A,0),(A,0)e£ + }, (12) 

where J is given by 

J(A, 0) = tr (AT n ) - trlog (E n kE T n + XJ N QXJ T N - 0) . (13) 

It can be shown (flU Theorem 6.1]) that the function J admits a unique minimum point in (A, 0) 
in C + . 

In this section we propose a modified gradient descent algorithm with backtracking line search 
(see, e.g., Ch. 9]) for the numerical solution of the dual problem ( [T2] ). This task requires some 
care because we are working in a matricial space. Let Yi^ N denote the orthogonal projection 
onto the linear subspace of symmetric, block-circulant matrices Cat. Before proceeding, we need 
two preliminary lemmas. The first one, gives a close form expression for the projection of 
E n AE^ onto <£n, while the second provides a useful characterization of matrices belonging to 
the orthogonal complement of €n- We refer the reader to flU for the proof of these statements. 

Lemma 4.1: Let A e & n +\ be the matrix 



A 



Ann A 



01 



A<ji A n 



A T A T 



Ann 

Ai„ 
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The orthogonal projection of E n AE^ onto £at, say IIa, is given by 



n A := U €n {E n KEl) 



n 


n7 


nj 




ni 


ni 


n 


n7 




n 2 


nj 




ni 


n 


n7 



n^ n^ ... n x n c 



with 



n = — (Aoo + An + • • • + Ar 



III 



n, 



A' 



(A 



01 



A 



12 



A ^l T 



1a t 

JV 0n ' 



(14a) 
(14b) 
(14c) 
(14d) 



while Hi = 0, forall i in the interval n + l<i<N — n — 1 . 

Lemma 4.2: Lef M <G ©at. 77ien M e (£n) ± if and only if it can be expressed as 



M = \J N N\Jl - N 



'N 



(15) 



for some N e &n- 

Consider now the functional 



(16) 



T 

J -n ■ 



J(A) := tr (AT n ) - tr log {n £jv (^A^) } 

whose gradient Va</(A) is given by 

V A J(A) = -El [U Cn {E n kEl)Y l E n 

The proposed algorithm is as follows. 

Theorem 4.1: Algorithm^is a gradient descent algorithm restricted to the subspace 

{(A, 9) | II,x (E n AEj t ) = - (U N QUj, - 0) } . 

Proof: Let (A, 6) be the unique minimum point of the functional J on C + . We know that 
(A, 0) are such that S° = E n AEj + U^QUj} — 6 is circulant. Thus, one can think of restricting 
the search for the solution of the optimization problem to the set 

{(A, 0) | (E n \El + U N BUjf - 0) is circulant} 



(17) 
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Algorithm 3 Matricial gradient descent algorithm 



Given a starting point A E dom J, a E (0,0.5), /3 E (0, 1) 
while ||V A J(A)|| 2 >r/do 

AA := -VaJ(A) 

while J(A + iAA) > J(A) + attr {VJ(A) T AA} do 
t :=/3t 

end while 

A := A + £AA 
end while 



which can be represented as 

{(A, 6) | n € x (E n kEl + U N QUl - 0) = o} 



Since {U N QUjj — 6) E <tjj (see Lemma 4.2), this expression can be rewritten as 



{(A, 9) | n € x (E n kEl) = - (UsQUl - e) } . 



If we compute the dual function J on the set ( fT7[ ), we obtain 

j(A,e) 



(A,e)|n e j.(s„ABT) = _( [/jv0[/ T_ e )J 

= tr (AT n ) - tr log (^A^ + ^0^ - T ) 

= tr (AT,) - tr log (E n AEj - n € x (E n AE^ 

= tr (AT n ) - trlog (n £jv (S„A^)) (18) 

which is the modified functional defined above. This shows that the proposed algorithm is nothing 
but a gradient descent algorithm in which the search of the minimum point has been restricted 



to the subspace where the optimal solution is known to lie, i.e. subspace (17). 



A. Initialization 

The following proposition provides a good starting point for the iterative procedure of Algo- 
rithm [31 
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Proposition 4.1: Let the block-To eplitz matrix made of the first n + 1, mxm covariance lags 

{So, Si, . . . , E n } 

T n := Toepl (S , Ei, . . . , S n ) , 

be positive definite and let {E^, k = 0, 1, 2, . . .} with 

E fc = E fc , fc = 0,1,2, ...,n 

be the maximum entropy (positive) extension of {E , Ei, . . . , E n }. Then, for N large enough the 

(c) 

block-circulant matrix Ey given by 
Toepl (E , E7,...,E^,EL 1 ,...,E^_ ,E^ + L, Ejv_ 1; . . . ,£ n +i, E n , . . . Ei ) , for N even 

\ 2 L 2 2 2 / 

Toepl (e , E7, • • • , E^", E^ +1 , . . . , E^-i , Ejv-l, . . . , E n+1 , E n , . . . , Y,A , for N odd 

(19) 

is a covariance matrix which for N — > oo is arbitrarily close to the mN x mN maximum 

entropy block-circulant extension of T n . 

Proof: That Ey is a valid covariance matrix for N large enough follows from [4, Theorem 



<(c) 



5.1]. To show that Xy given by ( |T9| ) tends to the maximum entropy block-circulant completion 



of the given data, it suffices to show that its inverse tends to be banded block-circulant for 
N — > oo (i.e. its off-diagonal 
can be block-diagonalized as 



(c) 

N — > oo (i.e. its off-diagonal blocks tends uniformly to zero). To this purpose, recall that Ey 



E$ = V*jvV* (20) 

where V is the Fourier block-matrix whose k, /-th block is 

V kl := 1/v^Nexp [-j2ir{k - 1)(? - 1)/N] I m (21) 

and \Ev is the block-diagonal matrix 

*Ar:=diag(* ,*i,---,*JV-i), (22) 

whose diagonal blocks ty e , are the coefficients of the finite Fourier transform of the first block 
row of Yj$ 

^ £ = E + <P'±1 + (e^) 2 Ej + • • • + (e^)^ 2 E 2 + (e^)" -1 E 1; (23) 

with $£ := —2tc£/N. Thus in particular 
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where 

*^:=diag(V,*iV--^-i)- 

Now, let $(z) := £ + X^i ^i z ~ % + ( SSi ^* z ~* ) be the spectral density matrix of the 
maximum entropy completion of T n . It is well-known ll30ll that &(z) can be expressed in factored 
form as 

*(z) = [L^z- 1 )]- 1 K [L n {z- l )Y* (24) 

where L n (z~ l ) is the ra-th Levinson-Whittle matrix polynomial associated with the block- 
Toeplitz matrix T n 

n 

L{z- l ) = Y,Mk)z- k (25) 

k=0 

with the A n (k)'s and A n = A^ > being the solutions of the Yule- Walker type equation 



A n (0) A n (l) . . . A n {n)\ T; = [A n ... 
Note that $(z) _1 = L n (z~ 1 )* A" 1 L n (z) is a Laurent polynomial, that can be written as 
^(z)' 1 = M + (M x z + M 2 z 2 + ■■■ + M n z n ) + (Mi* + M 2 z 2 + ■ ■ ■ + M n z n )* 
Moreover, since e~^ eN = 1, W, the \EVs can be written as [j 

^ £ = So + e^±J + ■■■+ {e^) h ± T h + e~ j ^Ei + • • • + (e~^) h t h 
where 



7V-1 



., . N odd 
N/2, N even 



Hence 



where 



ty £ = $ ( e J^) - [A$ w (e j ^) + A$^ (e j ^)] 



A$ N (z 



£* 



=h+l 



Since the causal part of $(z) is a rational function with poles inside the unit circle, 

sup \\A$ N (e j ^) + A^ (e j ^) II ->• 

i=0,...,JV-l 



(26) 



(27) 



(28) 



2 Notice that for AT even e^ = e"^^ = -1, so that (e^)' 1 tj + (e"^)' 1 S h = - (th + tj 
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exponentially fast for N —$■ oo. It follows that for N — > oo 

*7 X ->■ ($(e^)) _1 = M + M ie j ^ + • • • + M n (e j ^) n + M^ (^) N ~ 1 + ■ ■ ■ + Mj (e j ^) 



N-n 



(29) 



for £ = 0, 1, . . . , N — 1, i.e. \I/ £ * tends to the finite Fourier transform of a sequence of the form 



i.e. £ 



.(c) 



M ,M 1 T ,M 2 T ,...,M r [,0,...,0,M n ,...,M 1 , 



tends to be banded block-circulant, as claimed. 



An alternative way to compute the maximum entropy completion of a partially specified 
block-Toeplitz matrix is via the formula [[Toll . 



-l _ 



where 



with 



*(z) = G*(z)T~ l B (B"r?B) B*T~ l G(z 



G{z) 



zl -A) B 



B 



A 






I 











I 

















1 1 

It follows that the spectral factor W(z) := [L n (z~ 1 )} A« has a realization 

W(z) = C(zI -A)~ 1 B + D 



with D = AL C = - 



A 



— 


A n (n) A n (n - 


-1) • 




• A n (l) 





lm 

















-*m 









and 




-A n {n) 



-A n (n - 1) 



-A n {\) 



B 








A 2 

Jin 



(30) 



(3D 



(32) 
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The positive real part of the maximum entropy spectrum is given by 

$ + (z) = C(zI-A)- 1 C T + ^S (33) 

where C T = APC T + BD T , with P = APA T + BB T and the maximum entropy covariance 
extension results 

t k = CA k - l C T , k > n. 

With this extension at hand, we can compute an approximation for the maximum entropy block- 
circulant extension as suggested by Proposition |4T[ A good starting point for our gradient descent 
algorithm can then be obtained from ([14]) assuming for A a Toeplitz structure. 

V. Numerical experiments 

The matricial gradient descent algorithm has been implemented in Matlab. The results are 
shown in Figures [4] and [5] along with Tables [XT] and III The implementation exploits the block- 



circulant symmetric structure (recall, in particular, that for block-circulants the inverse can be 
computed efficiently by means of a Fourier transform). At each iteration, the algorithm requires 
the inversion of \^j^~\ matrices of order m. It follows that the execution time increases as the 
completion size N and the block size m increase (see Figure [4] and Table [TT]) . Finally, it also 
increases, although to a lesser amount, for larger bandwidth n (see Figure [5] and Table III). 



Table IV presents a comparison between the execution times for different starting point A . The 
gradient descent algorithm has been initialized to the normalized identity and by the procedure 
of Section IV-A| The proposed initialization acts effectively to reduce the number of iterations 



(and thus the computational time) to reach the minimum. 

Finally, the gradient descent algorithm (GD) has been compared to the iterative proportional 
scaling procedure (IPS) by Speed and Kiiveri. Both algorithms are implemented in Matlab. The 
Bron-Kerbosch algorithm [3 J has been employed for finding the cliques in the graph for IPS. The 
execution times for different completion size N and block size m are plotted in Figures [7] and [8] 
along with Tables |Vj and fVij It can be seen that the gradient descent algorithm runs faster than the 
iterative proportional scaling and the gap between the two increases as N increases. Moreover, 
the gap becomes much more evident as m grows, making the gradient descent algorithm more 
attractive for applications where the process under observation is vector-valued (m > 1). 
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Gradient descent algorithm - CPU time 




200 250 

N (completion size) 



Fig. 4: Matricial gradient descent algorithm: CPU time [sec.] for bandwidth n — 1, m — {1, 3}, 
and completion size N varying from 50 to 400. 



50 


0.5535 


0.9749 


100 


1.9376 


3.4989 


150 


4.2258 


7.7427 


200 

N 

250 


7.3857 
11.4440 


13.5903 
20.9953 


300 


16.3449 


30.0519 


350 


22.1412 


40.6536 



400 28.7854 52.7949 



TABLE II: Matricial gradient descent algorithm: CPU time [sec] plotted in Figure |4| for 
bandwidth n = 1, m = {1,3}, and completion size iV varying from 50 to 400. 
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Gradient descent algorithm - CPU time 




Fig. 5: Matricial gradient descent algorithm: CPU time [in sec.] for N = 50, m = 1, n varying 
between 2 and 20. 



n 


CPU time [sec] 


2 


1.1351 


4 


1.7895 


6 


1.9818 


8 


2.2215 


10 


3.4312 


12 


4.8058 


14 


5.2528 


16 


5.6626 


18 


7.3284 


20 


7.4922 



TABLE III: Matricial gradient descent algorithm: CPU time [in sec] plotted in Figure |5| for 
N = 50, m = 1, n varying between 2 and 20. 
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N 
10 

20 
30 
40 
50 



m 

5 
5 
5 
5 

5 



GDI 
# of itz. CPU time 



99 
212 
322 
432 
541 



0.1455 
0.4143 
0.8355 
1.4233 
2.1937 



GDT 
# of itz. CPU time 



61 

65 

97 

130 

163 



0.0767 
0.1270 
0.2504 
0.4285 
0.6603 



TABLE IV: CPU time [in sec] for the matricial gradient descent algorithm with different 



initializations (identity on the left and as in Section IV-A on the right). The reported times 
have been computed for N = [10, 20, 30,40, 50], m = 5 and bandwidth n — 3. 



Gradient descent - CPU time 




150 
Nm (ccmpletion size) 



Fig. 6: CPU time [in sec] for the matricial gradient descent algorithm with different initializations 



(identity in green and as in Section IV-A in blue). The reported times have been computed for 

N = [10, 20, 30, 40, 50], m = 5 and bandwidth n = 3. 
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TV 


m 


IPS 


GD 


10 


5 


4.7048 


0.0767 


20 


5 


16.4981 


0.1270 


30 


5 


29.2779 


0.2504 


40 


5 


43.8072 


0.4285 


50 


5 


63.8069 


0.6603 



TABLE V: Matricial gradient descent algorithm vs. iterative proportional scaling: CPU time [in 
sec] for N = [10, 20, 30,40, 50], m = 5, bandwidth n = 3. 





Gradient descent vs. Iterative proportional scaling - 


-CPU time 
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-GD, 




- 




* 


- IPS ?i 
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60 
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50 
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- 


£ 40 
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- 


CD 
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§30 


JT 






- 


20 


/^ 






- 


10 


/^ 


■jf 




J| 


5 


6 % % 

100 150 


200 




25 



Nm (completion size) 



Fig. 7: Matricial gradient descent algorithm vs. iterative proportional scaling: CPU time [in sec] 
for iV = [10, 20, 30, 40, 50], m = 5, bandwidth n = 3. 
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N 



IPS 



GD 



10 


10 


69.7671 


0.1516 


20 


10 


307.9596 


0.4459 


30 


10 


597.3791 


0.8988 


40 


10 


924.6431 


1.4798 


50 


10 


1341.0976 


2.2052 



TABLE VI: Matricial gradient descent algorithm vs. iterative proportional scaling: CPU time [in 
sec] for N = [10, 20, 30,40, 50], m = 10, bandwidth n = 3. 



Gradient descent vs. Iterative proportional scaling - CPU time 



1400 




100 



150 



200 



250 300 350 

Nm (completion size) 



400 



450 



500 



Fig. 8: Matricial gradient descent algorithm vs. iterative proportional scaling: CPU time [in sec] 
for N = [10, 20, 30, 40, 50], m = 10, bandwidth n = 3. 



VI. A NECESSARY AND SUFFICIENT CONDITION FOR THE FEASIBILITY OF THE CMEP FOR 

UNITARY BLOCK-SIZE AND BANDWIDTH ONE 

In this section, we provide a necessary and sufficient condition for the feasibility of the CMEP 
[4] for the special case of unitary block-size and bandwidth. We start by considering a real- valued, 
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discrete-time, stationary periodic process with circulant covariance matrix. We show that the 
(stochastic) DFT representation features uncorrelated elements. We then write explicitly the 
covariance lags of such a process. Next, we reformulate the feasibility of the CMEP for such 
a process in terms of the solvability of a system of linear equations and provide necessary and 
sufficient conditions for the solvability of this system. 

Theorem 6.1: Let {y(t)} be a zero mean stationary periodic process of period iV taking 
values in R. Then 

(i) {y(t)} can be represented as 

N-l 



where 

Cfc 



fc=0 

V-l 

N 



N-l 

^I>(*)e-^ 05) 



t=o 
and the c fc 's are uncorrelated random variables. 



(ii) the covariance samples of {y(t)} are given by 

^{e [|c | 2 ] + ^E [|c fc | 2 ]2cos(fcm 



K_ 1 

' 2n 



fc=i 



-E 



Cn 

2 



cos (rmr) > (36) 



for N even, and 



N-l 

^{E [|c | 2 ] +J2E [\c k \ 2 ] 2 cos (km^Pi } (37) 



fe=i 



for N odd. 



Proof: (i) Equations ([34]), ( [35] ) are simply the Discrete-Time Fourier Series equations for 
y(t). We want to show that the c^'s are uncorrelated random variables. To this aim, note that 
((34]) can be written in matrix form as 

y = V*c, (38) 



where V denotes the Fourier matrix (21) with m = 1, and y and c are the column vectors 



obtained by stacking the y(k)'s and the c^'s for k = 0, . . . , iV — 1. Using ( [38] ) and recalling that 
V is unitary, the variance of the vector c can be written as 

E [cc*] = VE N V\ (39) 
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where £ N := Eyy T . As we recalled in Section IV- A every circulant matrix can be diagonalized 



by the Fourier matrix. Since Tin is circulant, it follows from p9\ ) that E [cc*] is diagonal, i.e. 
the Cfc's are uncorrected random variables. 

(ii) By the uncorrelatedness of the c^'s, it follows that the covariance lags 



m 



a m = E [y(t + m)y(t)*} = E 
. . , N — 1, can be written as 






k(t+m)^ 



JV-1 



c e e ja2 N 



Or, 



\ X> [m 2 ] 



Jkm- 



N 



k=0 



Taking into account that c k 



-N- 



_ k , we get equations ( |36 l and (37). 



Theorem 6.1 let us reformulate the feasibility of the CMEP for m 



n 



1 in terms of the 



solution of a linear system of equations. This is the content of Corollary 6.1 below. 
Corollary 6.1: Let 

1 <7i 
Vl 1 

be the matrix of the boundary data, where, without loss of generality, we have assumed a = 1, 
and let N be a positive integer, N > 3. The data matrix T 2 admits a positive definite (respectively, 
a positive semidefinite) circulant completion if and only if the systems 



y [Po + J2k=i 2 Pk+Pf 

(— — l \ 

Po + Efc=i 2 Pfc cos (^#) ~Pf J = oi 



for N even 



(40) 



and 



iV-l \ 

F(Po + Efci 1 2p fc cos(A;f))=o- 1 



for N odd 



(41) 



with p k > (respectively, p fc > 0), k — 0, . . . , N — 1, admit solution. 



Proof: System ( |40| ) (respectively, system (|4T[)) follows simply by writing (|36j (respectively, 
( |3~7| )) for m = and m = 1. ■ 

We are now ready to state our main result. 
Theorem 6.2 (Necessary and sufficient condition - unitary block-size and bandwidth): 

The data matrix T 2 admits a positive definite circulant completion if and only if 
• [<7i| < 1, for N even; 
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cos (^tt) < 0-1 < 1, for N odd. 



Proof: Once again, by Corollary 6.1 it suffices to prove that system ( ]40[ ) for N even ((j4Tj) 
for iV odd) with the constraints pk > 0, k = 0, . . . ,N — 1 has solution if and only if |<7i| < 1 
(respectively, if and only if cos (^^^) < <J\ < 1)- For what concern necessity, notice that, for 
N even, since the p^s are > VA; 

-l 

/, 2tt\ 

Pn 



JV_i ^—1 

p q + 2_^2p k + pN >po + 2_^ 2 Pk cos I fc ^ J - 
fc=i fc=i ^ ' 



and indeed, since the p^'s must be strictly positive, the inequality is strict, i.e. 0\ < 1. Similarly, 

-l 

-Pn, 



(N-l \ ^-1 

Po+ 5Z 2pfc+p f <Po+ ^ 2 ^ cos ( fc jv ) 
fe=i / fe=i ^ ' 



so that ax > — 1 when the p^'s are strictly positive. Necessity for iV odd can be deduced by a 
similar argument. For what concern sufficiency, the two cases must be distinguish. 



Sufficiency for N even: we need to prove that if |cr x | < 1, then system (40) with constraints 

Pk > 0, k = 0, . . . , N — 1 has solution. In fact, setting 

p = N(a - e) 
Pn = JV(1 — a — s) 
with a G (0, 1), < e < min {a, 1 — a}, we get 

jf [po + Pk) = a — e + 1 — a — s = 1 — 2e 
jj [Po — Pk) = a — e — l + a + e = 2a — 1 
Thus if one chooses the remaining pk in order to satisfy 

(42) 
^E fc 2 = 1 1 2R.cos(A;f)=0 

— — l 
then (Jo = 1 and cti = 2a — 1 and the thesis is proved. Since ^ fc 2 =1 cos (^w) = ^' ^ * s eas y t0 

see that 




2 



-1 ' ' ' 2 



Pk — , N _, N ) & " " 1 5 • • • ; n " J 
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is a solution of d42l) and thus 



p = N(a - e) 
1 — a — e 

As 



Vn 

2 



Pk 



(f-D 



k 



N 
' 2 



-1 



with a = ^P", < e < a solves d40l) with constraints p k > 0, k — 0, . . . , N — 1. 



Sufficiency for iV odd: the aim is to prove that if cos (^j^^) < 0i < 1, then system (41 ) 

with constraints p k > 0, k = 0, . . . , iV — 1 has solution. If we set 
( 



p = N(a - e) + e 



P n-i 

2 



Pk 



Pe+i 

2 



N( 



S 



2 a ^2cos(V-) 

JV-1 jV+1 
' 2 ' 2 



+ 



A;^0 



+ e 



(43) 



with a G (0, 1), < e < min 



Na 



N 



a-1 



1-JV- 



' 2 



- + 1 



>, then the first of (41 ) is satis 



fled and the second yields <j\ — [l — cos (^p^)] ol + cos (^^7r), i.e. as a varies continuously 
over the interval (0, 1), <j\ varies continuously over (cos (^r^r) , l), which concludes the proof. 



Remark 6.1: Notice that a necessary and sufficient condition for the existence of a positive 
semidefinite circulant completion of the given data can be obtained by simply letting <j\ to take 
the values corresponding to the interval endpoints, i.e. the data matrix T 2 admits a positive 
semidefinite circulant completion if and only if |<Ji| < 1, for iV even, and cos (^r^) < Ci < 1, 
for JV odd. 

Remark 6.2: In their work ([20]) Grone, Johnson, Sa and Wolkowicz show that every partial 
positive definite matrix (i.e. every matrix such that all the specified principal submatrices are 
positive definite) has a positive definite completion if and only if the graph associated with the 
specified entries is chordal ([20, Theorem 7]). However, as we have already noticed, the graph 
associated with the given entries for the CMEP is not chordal and thus one may wonder if the 



result of Theorem 6.2 contrasts with the one in [EO 



The distinction becomes clear if the necessary and sufficient condition in 11201 is rephrased as 
follows: every partial positive definite matrix has a positive definite completion for every fixed 
completion size N if and only if the graph of the specified entries is chordal. On the other side, 
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Theorem |6.2| claims that, even for the non-chordal graph associated with the index set X h ([TJ) 
a positive definite completion of the given data may exist provided that the completion size is 
sufficiently large. So there is no contradiction between the two statements. 

VII. Example 

Example 7.1: Let 

1 0.01 

-0.91 1 
We want to investigate the feasibility of Problem [4] for N = 7 and N = 9. Since 

(N - 1) 1 -0.9010 for N = 7 

} -0.9397 for N = 9 



by Theorem 6.2 we expect that for iV = 7 the problem is unfeasible while for N > 9 it is 



expected to become feasible. For N = 7 the set of all positive definite completions is delimited 
by the intersection of the halfplanes delimited by (flU) 

ty( w °) = -0.82 + 2x + 2y (44a) 

^(w 1 ) = #( w 6 ) = -0.134751 - 0.445042a; - 1.80194y (44b) 

^(V ) = ^( w 5 ) = 1.40499 - 1.80194a; + 1.2469% (44c) 

^/( w 3 ) = ^(V) = 2.63976 + 1.24698s - 0.445042y. (44d) 

In Figure [9] the intersection T of the half-planes identified by (j44aj) and ( |44b[ ) is shown, together 
with the half-plane delimited by ( |44c[ ).The intersection of these two regions is empty. It follows 
that the intersection of the four half-planes is also empty, as claimed. On the other hand, if 
N = 9, the eigenvalues are 

y( w °) = -0.82 + 2x + 2y + 2z 

^(V ) = ^( w 8 ) = -0.394201 + 0.347296x -y- 1.879392 

^(V ) = ^(V ) = 0.68396 - 1.87939a; -y+ 1.532092 

\$j( w 3 ) = ^(w 6 ) = 1.91 - x + 2y - z 

^(u, 4 ) = ^<y ) = 2.71024 + 1.53209a; -y + 0.347296z 



and the feasible set is the nonempty region shown in Figure 10 
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y 
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Fig. 9: One half -plane and the intersection of other two of the four half-planes representing the 
regions where the eigenvalues of Circ {1, —0.91, x, y, y, x, —0.91} are positive. 




Fig. 10: Feasible region {(x, y, z) \ *E N > 0} for S N = Circ {1, —0.91, x, y, z, z, y, x, —0.91} 
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VIII. Conclusions 

The main contribution of the present paper is an efficient algorithm to solve the maximum 
entropy band extension problem for block-circulant matrices. This problem has many applica- 
tions in signal processing since it arises in connection with maximum likelihood estimation of 
periodic, and in particular quasi-Markov (or reciprocal), processes. Even if matrix completion 
problems have gained considerable attention in the past (think for example to the covariance 
extension problem for stationary processes on the integer line, i.e. for Toeplitz matrices), the 
maximum entropy band extension problem for block-circulant matrices has been addressed for 
the first time in [4]. The proposed algorithm heavily exploits the circulant structure and relies 
on the variational analysis brought forth in |3). Further light is shed also on the feasibility issue 
for such a problem. 
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